home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / dbesy.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  5.3 KB  |  148 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((nulim (make-array 2 :element-type 'f2cl-lib:integer4)))
  12.   (declare (type (simple-array f2cl-lib:integer4 (2)) nulim))
  13.   (f2cl-lib:fset (f2cl-lib:fref nulim (1) ((1 2))) 70)
  14.   (f2cl-lib:fset (f2cl-lib:fref nulim (2) ((1 2))) 100)
  15.   (defun dbesy (x fnu n y)
  16.     (declare (type double-float x fnu)
  17.              (type f2cl-lib:integer4 n)
  18.              (type (simple-array double-float (*)) y))
  19.     (prog ((w (make-array 2 :element-type 'double-float))
  20.            (wk (make-array 7 :element-type 'double-float)) (azn 0.0) (cn 0.0)
  21.            (dnu 0.0) (elim 0.0) (flgjy 0.0) (fn 0.0) (ran 0.0) (s 0.0) (s1 0.0)
  22.            (s2 0.0) (tm 0.0) (trx 0.0) (w2n 0.0) (xlim 0.0) (xxn 0.0) (i 0)
  23.            (iflw 0) (j 0) (nb 0) (nd 0) (nn 0) (nud 0))
  24.       (declare (type f2cl-lib:integer4 nud nn nd nb j iflw i)
  25.                (type (simple-array double-float (7)) wk)
  26.                (type (simple-array double-float (2)) w)
  27.                (type double-float xxn xlim w2n trx tm s2 s1 s ran fn flgjy elim
  28.                 dnu cn azn))
  29.       (setf nn (f2cl-lib:int-sub (f2cl-lib:i1mach 15)))
  30.       (setf elim (* 2.303 (- (* nn (f2cl-lib:d1mach 5)) 3.0)))
  31.       (setf xlim (* (f2cl-lib:d1mach 1) 1000.0))
  32.       (if (< fnu 0.0) (go label140))
  33.       (if (<= x 0.0) (go label150))
  34.       (if (< x xlim) (go label170))
  35.       (if (< n 1) (go label160))
  36.       (setf nd n)
  37.       (setf nud (f2cl-lib:int fnu))
  38.       (setf dnu (- fnu nud))
  39.       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 nd)))
  40.       (setf fn (- (+ fnu n) 1))
  41.       (if (< fn 2.0) (go label100))
  42.       (setf xxn (/ x fn))
  43.       (setf w2n (- 1.0 (* xxn xxn)))
  44.       (if (<= w2n 0.0) (go label10))
  45.       (setf ran (f2cl-lib:fsqrt w2n))
  46.       (setf azn (- (f2cl-lib:flog (/ (+ 1.0 ran) xxn)) ran))
  47.       (setf cn (* fn azn))
  48.       (if (> cn elim) (go label170))
  49.      label10
  50.       (if (< nud (f2cl-lib:fref nulim (nn) ((1 2)))) (go label20))
  51.       (setf flgjy -1.0)
  52.       (multiple-value-bind
  53.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
  54.           (dasyjy #'dyairy x fnu flgjy nn y wk iflw)
  55.         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6))
  56.         (setf iflw var-7))
  57.       (if (/= iflw 0) (go label170))
  58.       (if (= nn 1) (go end_label))
  59.       (setf trx (/ 2.0 x))
  60.       (setf tm (/ (+ fnu fnu 2.0) x))
  61.       (go label80)
  62.      label20
  63.       (if (/= dnu 0.0) (go label30))
  64.       (setf s1 (dbesy0 x))
  65.       (if (and (= nud 0) (= nd 1)) (go label70))
  66.       (setf s2 (dbesy1 x))
  67.       (go label40)
  68.      label30
  69.       (setf nb 2)
  70.       (if (and (= nud 0) (= nd 1)) (setf nb 1))
  71.       (dbsynu x dnu nb w)
  72.       (setf s1 (f2cl-lib:fref w (1) ((1 2))))
  73.       (if (= nb 1) (go label70))
  74.       (setf s2 (f2cl-lib:fref w (2) ((1 2))))
  75.      label40
  76.       (setf trx (/ 2.0 x))
  77.       (setf tm (/ (+ dnu dnu 2.0) x))
  78.       (if (= nd 1) (setf nud (f2cl-lib:int-sub nud 1)))
  79.       (if (> nud 0) (go label50))
  80.       (if (> nd 1) (go label70))
  81.       (setf s1 s2)
  82.       (go label70)
  83.      label50
  84.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  85.                     ((> i nud) nil)
  86.         (tagbody
  87.           (setf s s2)
  88.           (setf s2 (- (* tm s2) s1))
  89.           (setf s1 s)
  90.           (setf tm (+ tm trx))
  91.          label60))
  92.       (if (= nd 1) (setf s1 s2))
  93.      label70
  94.       (f2cl-lib:fset (f2cl-lib:fref y (1) ((1 *))) s1)
  95.       (if (= nd 1) (go end_label))
  96.       (f2cl-lib:fset (f2cl-lib:fref y (2) ((1 *))) s2)
  97.      label80
  98.       (if (= nd 2) (go end_label))
  99.       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
  100.                     ((> i nd) nil)
  101.         (tagbody
  102.           (f2cl-lib:fset (f2cl-lib:fref y (i) ((1 *)))
  103.                          (-
  104.                           (* tm
  105.                              (f2cl-lib:fref y
  106.                                             ((f2cl-lib:int-sub i 1))
  107.                                             ((1 *))))
  108.                           (f2cl-lib:fref y ((f2cl-lib:int-sub i 2)) ((1 *)))))
  109.           (setf tm (+ tm trx))
  110.          label90))
  111.       (go end_label)
  112.      label100
  113.       (if (<= fn 1.0) (go label110))
  114.       (if (> (* (- fn) (- (f2cl-lib:flog x) 0.693)) elim) (go label170))
  115.      label110
  116.       (if (= dnu 0.0) (go label120))
  117.       (dbsynu x fnu nd y)
  118.       (go end_label)
  119.      label120
  120.       (setf j nud)
  121.       (if (= j 1) (go label130))
  122.       (setf j (f2cl-lib:int-add j 1))
  123.       (f2cl-lib:fset (f2cl-lib:fref y (j) ((1 *))) (dbesy0 x))
  124.       (if (= nd 1) (go end_label))
  125.       (setf j (f2cl-lib:int-add j 1))
  126.      label130
  127.       (f2cl-lib:fset (f2cl-lib:fref y (j) ((1 *))) (dbesy1 x))
  128.       (if (= nd 1) (go end_label))
  129.       (setf trx (/ 2.0 x))
  130.       (setf tm trx)
  131.       (go label80)
  132.      label140
  133.       (xermsg "SLATEC" "DBESY" "ORDER, FNU, LESS THAN ZERO" 2 1)
  134.       (go end_label)
  135.      label150
  136.       (xermsg "SLATEC" "DBESY" "X LESS THAN OR EQUAL TO ZERO" 2 1)
  137.       (go end_label)
  138.      label160
  139.       (xermsg "SLATEC" "DBESY" "N LESS THAN ONE" 2 1)
  140.       (go end_label)
  141.      label170
  142.       (xermsg "SLATEC" "DBESY" "OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL" 6
  143.        1)
  144.       (go end_label)
  145.      end_label
  146.       (return (values nil nil nil nil)))))
  147.  
  148.